TWENTY DIGITS OF SOME INTEGRALS OF THE PRIME 

ZETA FUNCTION 



RICHARD J. MATHAR 



Abstract. The double sum Y^ S >1 J2 P 1 /(P° l°gP s ) = 2.00666645 . . . over the 
inverse of the product of prime powers p s and their logarithms, is computed 
to 24 decimal digits. The sum covers all primes p and all integer exponents 
s > 1. The calculational strategy is adopted from Cohen's work which basically 
looks at the fraction as the underivative of the Prime Zeta Function, and then 
evaluates the integral by numerical methods. 



1. Overview 

Definition 1. The constant in the focus of this work is the sum 



(!) C^VV-^- 

over all primitive prime powers p s of the prime numbers p \15\ 2.27.3] [T4l §2.2]. 

Remark 1. The prime powers are represented by Shane's sequence A000961 — we 
drop the leading term p s — 1 to avoid division by zero — . The constant C and the 
contribution of s = 1 are entry A137250 and A137245 [26] . 

The only aim of the work is to improve on the previous estimates 2.00 < C < 2.01 
[12] and C > 2.006 0]. 

The simple computational strategy is to accumulate the partial sums over s in 

oo 

( 2 ) c^Y,-Y.^—> 

^ s ^ p s logp 

s=l p ° 

which converge at reasonable speed to an accuracy of 10 -7 in C after 20 terms or 
to an accuracy of 10 -19 after 59 terms, for example. 



2. Cohen's Integral 

The implementation is a shameless replica of Henri Cohen's reduction of the 
double sum over exponents and primes to a series of integrals [S] . 
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2.1. Logarithm-to-Integral Conversion. The logarithm of Euler's formula for 
Riemann's zeta function [3] [HI 9.523.1] 



(3) f(«) = IIr^ 

- L - L 1 — n 



is 



(4) logC(«) = -J>g(l-p-), 

p 

which turns with the Taylor expansion of the logarithm into [T51 9.523.2] j 

fe=l p y 

Definition 2. T/ie Prime Zeta Function is IT6l l23l 



(6) p (*) = Ei Ks>1 



p 



So the penultimate equation can be rephrased as 

(7) iogC(s) = £r p (M- 

fc=i 

The Mobius inversion of this formula reads [H (9b)][ini §17.1.3] 

(8) P(s ) = g^)i ogC(fcs ). 

k=l 

First integrals of terms in ((6|) are c 9, 565.1.] 

1 , 1 
— dt = — — . 

p l p s log p 



(9) / ~ [ dt=- 7T — 1 s>l, 



so integration of ([5]) using © gives 
(10) 

Definition 3. Cohen's integral is 



(11) J(m) = / logC(*)d* 

«/ m 

over the logarithm of Riemann's zeta function with variable integer lower limit [5]. 
Insertion into (fTU|) outlines the strategy to evaluate ©, 

(1 2 ) E-r— = E# J ( Jfe *)- 

v ; ^p s logp ^ /c 2 

p r b ^ fc=l 
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2.2. Transition to the Dirichlet Eta Function. For 1(1), the pole of at 
t = 1 can be isolated by handling the singularity via the smooth function [TJ 
23.2. 19] p 267] 

(13) r,(a) = (l-2 1 -')C(s). 

For lower limits of the integral larger than one, this recipe is not needed; we stick 
to it to present a shorter, simpler program. 

Remark 2. The finite value at s = 1 is [Tg][18l 0.232.1,1.511] 

(14) rj(l) = log 2. 

This follows also from the residuum of Res |s=i = 1 in conjunction with the 

Taylor expansion 1-2 1 -* = - ££l ( ~ l °f 2) ' (s - 1)' = (s - 1) log 2 - ^ (s - l) 2 + 
0((s — l) 3 )- //we introduce Stieltjes constants jj [SI® [7], 

(15) C(s) = _l_ +7 _ 7l(s _i) + ^ (s _i)2 + ... ^ 

i/ie Taylor series of r\ near s = 1 becomes 
(16) 

*W = ^ 2 + log 2 ( 7 - ^ ) (s - 1) - log 2 ( 71 + 7^ - ^) (- - I) 2 + • • • • 
The integrated logarithm of JT5J) is 

pOO poo pOO 

(17) / ds\ogn(s) = dslog(l-2 1 - s )+ / dslogC(s). 
One term is a dilogarithm [TTl [T^l |2"U] . 

(18) ^ \og(l-2 1 - s )ds= [ - l °f X dx = --^—U 2 



log2(l-x) log 2 

Remark 3. Special values at m = 1 and m — 2 are Li2(l) = 7r 2 /6 and Li2(l/2) = 
tt 2 /12 - ±log 2 2; see the constants A013661 and A076788 in the On-Line Ency- 
clopedia of Integer Sequences [26]. An accurate representation of ir 2 / (12 log 2) is 
entry A100199. 77(1) = log 2 is A002162, 7 is A001620, r( (1) = log 2(7 - log 2/2) 
zs A091812, and -71 is A0826S3. 

This moulds US into 



(19) J(m) = / ds log 77(5) + — Li 2 



>m— 1 



2.3. Numerical Implementation. The variable substitution 

m m m 

(20) u = l ; s—- ; ds — — rrrdu 

s 1 — u (1 — u) z 

maps the interval m < s < 00 onto the interval < u < 1. 

Remark 4. Alternative substitutions like u — 1/m — m/s 2 or u = 1 — 2~ s also 
map the half-infinite s-interval to a finite u-interval. They have the advantage of 
more evenly balanced integral kernels, but the disadvantage of infinite slopes at one 
end of the u-interval. 
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Figure 1 . The functions (m log rj)/(l - u) 2 < or (m log 0/(1 - 
u) 2 > arising in the integral kernels through the substitution 

(HOD- 



Wynn's e-algorithm 28, 22 is applied to numerical values gathered by the trape- 
zoidal rule to evaluate the 77-integral in (fTTH) . 

r°° r 1 lo ^v (t^i) 

(21) / ds\ogr)(s) =m / v ' du. 

Jm JO I 1 — u ) 

Remark 5. Because logC(s) « 2 _s at s > 10, I(m) w l/(2 m log2) as m — » 00. 
T/iis functional dependence leads to an almost straight line on the semi-logarithmic 
plot in Fig. [H So the first neglected term in (0) is a good estimator to the error in 
the partial sums. 

2.4. Intermediate Results. The following table shows s and I(s) in the left dou- 
ble column, s and ^2 p l/(p s logp) in the right double column, last digits rounded: 



1 


1 


797569958628739407930251e+00 


1 


1 


636616323351260868569658e+00 


2 


5 


365269459211771009617190e- 


-01 


2 


5 


077821878591993187743751e- 


-01 


3 


2 


275044530363196437654595e- 


-01 


3 


2 


212033403969814969599433e- 


-01 


4 


1 


041502253168599126451383e- 


-01 


4 


1 


026654782331571962664758e- 


-01 


5 


4 


942465209576037094607771e- 


-02 


5 


4 


906357474830439568004180e- 


-02 


6 


2 


392429939526549813915300e- 


-02 


6 


2 


383519825459824737569112e- 


-02 


7 


1 


171800374394507398697301e- 


-02 


7 


1 


169586557789515115400634e- 


-02 


8 


5 


781462903507255601662420e- 


-03 


8 


5 


775944595316604352784204e- 


-03 


9 


2 


865717415930253367688719e- 


-03 


9 


2 


864339771659922115926198e- 


-03 


10 


1 


424706500149304604154559e- 


-03 


10 


1 


424362320209387639441021e- 


-03 


11 


7 


096782687818267364308355e- 


-04 


11 


7 


095922515383712411613441e- 


-04 
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Figure 2. The functions I(s) and pS * fall of exponentially 
as a function of their argument s, and are not distinguishable in 
this common plot. 



12 


3 


539573278954885950994965e- 


-04 


12 


3 


539358269254206683875757e- 


-04 


13 


1 


766870488681893665136623e- 


-04 


13 


1 


766816740292220402778692e- 


-04 


14 


8 


824687187128450251937842e- 


-05 


14 


8 


824552821042831521626792e- 


-05 


15 


4 


409135113515804839684079e- 


-05 


15 


4 


409101522588177584140913e- 


-05 


16 


2 


203501275340515891846297e- 


-05 


16 


2 


203492877680910955374197e- 


-05 


17 


1 


101395992552518507579192e- 


-05 


17 


1 


101393893146441713313451e- 


-05 


18 


5 


505799725855358158700238e- 


-06 


18 


5 


505794477350959773416061e- 


-06 


19 


2 


752506920732716039062161e- 


-06 


19 


2 


752505608607939317868454e- 


-06 


20 


1 


376122595515799546777009e- 


-06 


20 


1 


376122267484767801052807e- 


-06 


21 


6 


880177047831559600925881e- 


-07 


21 


6 


880176227754180030145763e- 


-07 


22 


3 


439943284950306802000172e- 


-07 


22 


3 


439943079930986522426868e- 


-07 


23 


1 


719923247093458602614288e- 


-07 


23 


1 


719923195838631569242761e- 


-07 


24 


8 


599454961284003932811501e- 


-08 


24 


8 


599454833146940100459182e- 


-08 


25 


4 


299673733467609442174442e- 


-08 


25 


4 


299673701433343948018079e- 


-08 


26 


2 


149818953720244919068455e- 


-08 


26 


2 


149818945711678602970514e- 


-08 


27 


1 


074903506531325106214723e- 


-08 


27 


1 


074903504529183534309235e- 


-08 


28 


5 


374497633245768258586266e- 


-09 


28 


5 


374497628240414337653354e- 


-09 


29 


2 


687242183906380835200363e- 


-09 


29 


2 


687242182655042356063437e- 


-09 


30 


1 


343618881152634756065114e- 


-09 


30 


1 


343618880839800136417081e- 


-09 


31 


6 


718087036690018174171263e- 


-10 


31 


6 


718087035907931625220498e- 


-10 


32 


3 


359041062052731307704559e- 


-10 


32 


3 


359041061857209670487929e- 


-10 


33 


1 


679519712278619904984955e- 


-10 


33 


1 


679519712229739495683419e- 


-10 


34 


8 


397595832274787450826827e- 


-11 


34 


8 


397595832152586427576250e- 


-11 


35 


4 


198797006441494468678066e- 


-11 


35 


4 


198797006410944212865828e- 


-11 
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36 2 . 099398199991330939924248e-ll 

37 1.049698998919830689452100e-ll 

38 5 . 248494657681297045976310e-12 

39 2 . 624247216535094520600355e-12 

40 1 . 312123570832462108601650e-12 

41 6 . 560617729378942125931467e-13 

42 3 . 280308823095077105921359e-13 

43 1.640154397682756111408190e-13 

44 8.200771942197877942789073e-14 

45 4. 100385955693647813828150e-14 

46 2 . 050192972711729282745366e-14 

47 1 . 025096484644167040382940e-14 

48 5 . 125482417515178049571403e-15 

49 2 . 562741206855703686734819e-15 

50 1 . 281370602793890158864587e-15 

51 6.406853011856245416428711e-16 

52 3 . 203426505223720974861043e-16 

53 1 . 601713252377059924465468e-16 

54 8.008566261102631116158675e-17 

55 4 . 004283130290426065285105e-17 

56 2 . 002141565058249870693295e-17 

57 1 . 001070782500137215275764e-17 

58 5 . 005353912404060344256377e-18 

59 2 . 502676956169821595115832e-18 

60 1 . 251338478074174605310912e-18 
616. 256692390335085719290662e-19 

62 3 . 128346195155613757280562e-19 

63 1 . 564173097573830511199490e-19 

64 7 . 820865487855897997896814e-20 

65 3 . 910432743923530812923694e-20 

66 1 . 955216371960292677789151e-20 

67 9 . 776081859796554293375620e-21 

68 4 . 888040929896640781499146e-21 

69 2 . 444020464947774935687030e-21 

70 1 . 222010232473705649489420e-21 

71 6 . 110051162367922186267000e-22 

72 3 . 055025581183759072740187e-22 

73 1 . 527512790591812196239003e-22 

74 7 . 637563952958836514091412e-23 

75 3 . 818781976479343434677847e-23 

76 1 . 909390988239646776549639e-23 

77 9 . 546954941198150746783920e-24 

78 4 . 773477470599047661403870e-24 

79 2 . 386738735299514593372572e-24 

80 1 . 193369367649754217576498e-24 



36 2 . 099398199983693375971239e-ll 

37 1.049698998917921298463855e-ll 

38 5 . 248494657676523568505703e-12 

39 2 . 624247216533901151232704e-12 

40 1 . 312123570832163766259737e-12 

41 6 . 560617729378196270076685e-13 

42 3 . 280308823094890641957663e-13 

43 1 . 640154397682709495417266e-13 

44 8.200771942197761402811764e-14 

45 4. 100385955693618678833822e-14 

46 2 . 050192972711721998996784e-14 

47 1 . 025096484644165219445795e-14 

48 5 . 125482417515173497228539e-15 

49 2 . 562741206855702548649103e-15 

50 1 . 281370602793889874343158e-15 

51 6.406853011856244705125138e-16 

52 3 . 203426505223720797035150e-16 

53 1 . 601713252377059880008995e-16 

54 8.008566261102631005017492e-17 

55 4 . 004283130290426037499809e-17 

56 2 . 002141565058249863746971e-17 

57 1 . 001070782500137213539183e-17 

58 5 . 005353912404060339914924e-18 

59 2 . 502676956169821594030469e-18 

60 1 . 251338478074174605039571e-18 
616. 256692390335085718612310e-19 

62 3 . 128346195155613757110974e-19 

63 1 . 564173097573830511157093e-19 

64 7 . 820865487855897997790822e-20 

65 3 . 910432743923530812897196e-20 

66 1 . 955216371960292677782526e-20 

67 9 . 776081859796554293359059e-21 

68 4 . 888040929896640781495005e-21 

69 2 . 444020464947774935685995e-21 

70 1 . 222010232473705649489161e-21 

71 6 . 110051162367922186266353e-22 

72 3 . 055025581183759072740025e-22 

73 1 . 527512790591812196238962e-22 

74 7 . 637563952958836514091310e-23 

75 3 . 818781976479343434677822e-23 

76 1 . 909390988239646776549633e-23 

77 9 . 546954941 198150746783904e-24 

78 4 . 773477470599047661403866e-24 

79 2 . 386738735299514593372571e-24 

80 1 . 193369367649754217576498e-24 



At large s, the two values equalize because the first term at k = 1 dominates the 
sum at the right hand side of (fT!?)) . Cohen reported the first two rows |S]. Erdos 
and Zhang published an upper bound of ^/(p^°&p) [H]i a proof of convergence 
had been given earlier [15] , 
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The rightmost column P L (s, 1) is an aid to computation of the more general 



(22) P L ( S ,a)^J2 



(a - 1 + p) s logp 



E 

i=o 



s + l-1 

s - 1 



The /-series converges for a = and a 
with the alternating geometric series 
(23) 

P L (s, o+l) = £ -i 

P (a-l+p) s (H 



a)'P L (Z + s,l)- 
2. Larger a are then reached recursively 



i+P 



logp 



OO 

E 

i=0 



( a ^ 1 1 )(-l)'p L (.+l,«). 



This yields the following short table with double columns of s, a and P L (s, a) each: 



1 





2 


. 564343220686309193e+00 


4 


2 


2 


, 200937333087681441e- 


-02 


2 





1 


. 735535173734295407e+00 


5 


2 


6 


. 924273989916642133e- 


-03 


3 





1 


. 569430040669505312e+00 


6 


2 


2 


. 216719348246691407e- 


-03 


4 





1 


. 502480490856796843e+00 


7 


2 


7 


. 177059901294852453e- 


-04 


5 





1 


. 471819234100400691e+00 


8 


2 


2 


, 341804937450834077e- 


-04 


6 





1 


. 457080819283874067e+00 


9 


2 


7 


. 683437427098209717e- 


-05 


7 





1 


. 449846098516752367e+00 


10 


2 


2 


.531100154045370884e- 


-05 


8 





1 


. 446260454845751749e+00 


11 


2 


8 


. 362846790350413196e- 


-06 


9 





1 


. 444475273576531254e+00 


12 


2 


2 


. 769232699593665178e- 


-06 


10 





1 


. 443584547482253809e+00 


13 


2 


9 


, 185072149314254133e- 


-07 


11 





i. 


,443139643195342866e+00 


14 


2 


3 


, 050306451883795961e- 


-07 


12 





l 


, 442917304533643634e+00 


15 


2 


1 


, 013929603798279887e- 


-07 


13 





l 


,442806163373811902e+00 


16 


2 


3 


. 372678630308857437e- 


-08 


14 





l 


. 442750599803600762e+00 


17 


2 


1 


, 122456404716334707e- 


-08 


15 





l 


. 442722819765431096e+00 


18 


2 


3 


, 737099942191877583e- 


-09 


16 





i. 


,442708930182166938e+00 


19 


2 


1 


. 244595144792613520e- 


-09 


17 





l 


,442701985499337981e+00 


20 


2 


4. 


. 145889249942980892e- 


-10 


18 





l 


, 442698513185098959e+00 


1 


3 


1 


, 083884561790478549e+00 


19 





l 


, 442696777034769093e+00 


2 


3 


1 


. 539997501094314544e- 


-01 


20 





l 


. 442695908961300868e+00 


3 


3 


3 


. 279427700685188571e- 


-02 


1 


2 


l 


. 280433764453964036e+00 


4 


3 


7 


. 457583410785408360e- 


-03 


2 


2 


2 


. 518771844123868571e-01 


5 


3 


1 


, 747751882461897344e- 


-03 


3 


2 


7 


, 208846839108258174e-02 













In the left column of this table we observe lim s ^oo P L (s, 0) = 1/ log(2) w 1.44 . . .. 

3. Summary 

We obtained the constant 
(24) C = 2.00666645283106875643229 . . . 

by a rather basic numerical approach to Cohen's integral over the logarithm of the 
Prime Zeta Function. A table of J2 P l/(p s l°gp) was generated for 1 < s < 80 and 
a table of J2 P V[(P ± X ) S lo grf for 1 < s < 20. 

Appendix A. PARI Program 

The following PARI/GP program [M] implements the algorithm. WynnEpsItr 
and WynnEps calculate Wynn's generalized Richardson extrapolation of a sequence. 
logEta computes log?7. ICohen calculates (|T9| . Ieta calculates the integral (j2"Tj) . 
logSum returns the sum (|12p for a given argument s. ErdosLConst returns a value 
of©. 
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/** Wynn's epsilon process. One step of the recurrence. 

* @param epsl A column in the scheme of deltas. 

* @param eps2 The column to the right of epsl. 

* (Dreturn The extrapolation of the column epsl. 
*/ 

WynnEpsItr (epsl , eps2) ={ 

/* eps3 is the new column to the right of eps2. 

* slen is the number of items in epsl. 
*/ 

local (eps3, slen) ; 
slen=length(epsl) ; 
eps3=vector(slen-2) ; 

/* one-by-one filling of the new column */ 
for(i=l,slen-2, 

eps3[i]= epsl[i+l]+l/(eps2[i+l]-eps2[i]) ; 

) ; 

/* Either we have reached the rightmost column, indicated by the 

* fact that only one entry is left in eps3, or we iterate 

* once more with information now in the two 

* rightmost columns, eps2 and eps3. 
*/ 

if ( length (eps3) ==1 , 
return ( eps3[l]) , 
return ( WynnEpsItr (eps2,eps3)) ; 

) ; 

> 

/** Wynn's epsilon process. Main entry. 

* @param S The vector of the current estimates with an odd number of terms. 

* @return The extrapolation of the vector terms to infinity. 
*/ 

WynnEps(S)={ 

/* Essentially build the 2nd column of the scheme from the 

* inverted first differences of the vector, and enter the iteration. 
*/ 

local(epsl,slen = length(S)) ; 
epsl=vector (slen-1) ; 
f or (i=l , slen-1 , 

epsl[i]=l/(S[i+l]-S[i]) ; 

) ; 

return ( WynnEpsItr (S , epsl) ) ; 

> 

/** Logarithm of eta(s) . 

* The eta function is the product of 1-2" (1-s) with the Riemann zeta 

* function zeta(s) . The routine is meant to deal explicitly with 

* the value eta(l)=log(2) at the pole of zeta(s) . 
*/ 

logEta(s)={ 
if (s==l, 

return ( log(log(2)) ) , 

return( log( (1-2" (1-s) ) *zeta(s) ) ) ; 

) ; 
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} 

/** Calculate Cohen's I(m) = int_ (s=m. . infinity) log zeta(s) ds. 

* @param m The lower limit of the integral 

* @param relerr The relative error admissible by the result . 
*/ 

ICohen(m , relerr ) ={ 
local(li2,iz) ; 

/* Calculate a value of the Li2(l/2~ (m-1) ) into li2. 
*/ 

if (m==l , 

li2=Pi~2/6 , 
if (m==2 , 

li2=Pi"2/12-(log(2))"2/2 , 

li2=dilog(l/2~(m-l)) ; 

) ; 

) ; 

/* Delegate the main calculation to the integral of the 

* eta-function. 
*/ 

iz = Ieta(m, relerr) + li2/log(2) ; 
printpO'I ",m," ",iz) ; 
return ( iz ) ; 

} 

/** Calculate the integral int_ (s=m. . infinity) log eta(s) ds . 

* Implemented as m * int_0~l log eta(m/(l-u)) du/(l-u)~2. 

* @param m The lower limit of the integral 

* @param relerr Admissible relative error of the result 

* (Dreturn The value . 
*/ 

let a (m, relerr )={ 

/* simpO is the value of the kernel at the lower limit u=0. u are specific 

* abscissa points. WynnL is the column number in Wynn's extrapolation, an odd 

* number equal to or larger than 1 . 
*/ 

local(simpO=logEta(m) ,u, sus ,N, otrap, WynnL=9, sim) ; 

/* sus[l,.] contains sums of the integral kernel of the previous 

* step size; sus [2,.] accumulates the sums at the doubled abscissa count N. 

* We keep track of the sums in two batches, sus[,l] and sus[,2] for 

* even and odd indexed abscissa points, because subdivision of the 

* old intervals in factors of 2 can use the old sums again. 
*/ 

sus=matrix(2,2) ; 

/* sim[l] contains the previous extrapolation to zero step size, 

* sim [2] the newer extrapolation estimated with N abscissa values. 
*/ 

sim=vector (2) ; 
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/* otrap: A history of evaluations of the integral with trapezoidal rules 

* (which ensures that Wynn's extrapolation is applicable because the error is linear 

* in the step width) . otrap [i+1] obtained at half the step width of otrap [i] . 
*/ 

otrap=vector (WynnL) ; 

N=4 ; 
while (1 , 

/* Accumulate sum over function values at new points (interstitial 

* to previous points) in su[2,2] . 
*/ 

sus[2,2]=0 ; 
f orstep(i=l ,N-1 ,2 , 
u = i/N ; 

sus[2,2] += logEta (m/(l-u)) /(l-u)~2 ; 

) ; 

/* Sum of function values over the other half of the abscissa points is 

* sum of values over all points of the previous coarser grid. 
*/ 

sus[2,l]=sus[l,l]+sus[l,2] ; 

/* If this is the first loop, su[2,l] is zero (from the missing 

* previous loop) and is computed explicitly. 
*/ 

if (sus[2,l]==0, 

forstep(i=2,N-2,2, 
u = i/N ; 

sus[2,l] += logEta(m/(l-u)) / (l-u)~2 ; 

) ; 

) ; 

/* Shift the history stack of old values with coarser abscissa by one index */ 
f or (h=2, WynnL, 

otrap [h-1] =otrap [h] 

) ; 

/* The integral is (trapez . rule) sum of function values, half weight 

* to values at both ends of the interval. Factor m from the variable substitution. 

* Factor 1/N represents the subinterval width. 
*/ 

otrap [WynnL] = m* (simpO/2+sus [2 ,2] +sus [2, 1] ) /N ; 
if ( otrap [1] !=0. , 

/* If history has reached sufficient depth defined by WynnL, compute 

* Wynn's extrapolation in sim[2] , and leave loop if converged 

*/ 

sim[2] = WynnEps (otrap) ; 

if( abs (sim [2] -sim[l] ) < relerr*abs(sim[2] ) , 
return(sim[2] ) ; 

) ; 

) ; 

/* Save sums over function values of both submeshes for next iteration. */ 
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sus[l,l]=sus[2,l] ; 
sus[l,2]=sus[2,2] ; 
sim[l] =sim [2] ; 

/* Double number of subintervals . Half abscissa mesh width. */ 
N *= 2 ; 

) ; 

} 

/** Compute sum_p l/(p~s*log(p) over the primes p. 

* @param s The exponent in the denominator 

* @param relerr The relative error of the result 

* @return The sum including all terms until the first one 

* skipped is smaller than the error bar. 
*/ 

logSum(s , relerr) ={ 

/* a and aold are the new and previous partial sum. 

* k is the summation index, thisl is the integral I(k*s) 

*/ 

local (a=0,aold=0,k=l, thisl) ; 
while (1 , 

/* To avoid premature termination of the summation 

* thru a zero value of the Moebius function, add a 

* check on the function in advance . 
*/ 

if( moebius (k) != 0, 

/* The value of I(k*s) is either in the hash 

* vector, or calculated and saved into the vector. 
*/ 

if( k*s > length ( Ihash) , 

thisI=ICohen(k*s, relerr) , 

if( Ihash [k*s]==0, 

Ihash [k*s] =ICohen(k*s , relerr) 

) ; 

thisl=lhash[k*s] ; 

) ; 

/* Update partial sum by adding mu(k)*I(k*s)/k~2 

* to the previous value . 
*/ 

a += moebius (k) *thisl/k~2 ; 

if( abs(a-aold) < relerr*abs (a) , 

printpC'log ",s," ",a) ; 

return (a) ; 

) ; 

) ; 

k++ ; 
aold=a ; 

) ; 

} 
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/** Calculate the main constant . 

* @param relerr The relative error admitted to the result . 

* Sums the individual terms of sum_{s=l ,2 ,3 . . . . } 1/s sum_-[p=2, 3, 5, 7, 11 . . . } l/(p~s log p) 
*/ 

ErdosLConst (relerr) ={ 

local (a=0,aold=0,s=l) ; 

/* A specific I() integral occurs multiple times as various values 

* of the summation variable s, so gathering the values computed in 

* a hash vector saves time. Ihash[s] stores I(s) with capacity for s=1..400. 
*/ 

Ihash=vector (400) ; 
while (1, 

/* we are adding terms until the last one at s=sL falls below 

* the error margin. The last one is of the order log( . . (sL) . ) /sL 

* where log(...) is roughly l/(2~sL log2) . We could allow a relative 

* error larger than the original call here, but skip this 

* aspect for simplicity, not to jeopardize the use of the 

* hashed table of I() integrals. 
*/ 

a += logSum(s, relerr) /s ; 

print ("partial sum ", s," ", a) ; 

/* The relative errors is roughly half the absolute error, 

* so we leave the loop if the last term is smaller than the error. 
*/ 

if( abs(a-aold) < relerr*abs (a) , 
return (a) ; 

) ; 

aold=a ; 
s++ ; 

) ; 

} 

/** Main program. */ 
{ 

ErdosLConst (l.e-28) ; 
quit() ; 

} 
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